APS preprint 



Power Law Distributions of Offspring and Generation Numbers in Branching Models 

of Earthquake Triggering 

A. Saichev, 1,2 A. Helmstetter, 2 and D. Sornette 2 ' 3 ' 4 

1 Mathematical Department, Nizhny Novgorod State University, 
Gagarin prosp. 23, Nizhny Novgorod, 603950, Russia 
^Institute of Geophysics and Planetary Physics, 
University of California, Los Angeles, CA 90095 
s Department of Earth and Space Sciences, University of California, Los Angeles, CA 90095 
* Laboratoire de Physique de la Matiere Condensee, 
CNRS UMR 6622 and Universite de Nice-Sophia Antipohs, 06108 Nice Cedex 2, Franc^\ 

(Dated: February 2, 2008) 

We consider a general stochastic branching process, which is relevant to earthquakes as well as 
to many other systems, and we study the distributions of the total number of offsprings (direct 
and indirect aftershocks in seismicity) and of the total number of generations before extinction. We 
apply our results to a branching model of triggered seismicity, the ETAS (epidemic-type aftershock 
sequence) model. The ETAS model assumes that each earthquake can trigger other earthquakes 
("aftershocks"). An aftershock sequence results in this model from the cascade of aftershocks of 
each past earthquake. Due to the large fluctuations of the number of aftershocks triggered directly 
by any earthquake ("fertility"), there is a large variability of the total number of aftershocks from 
one sequence to another, for the same mainshock magnitude. We study the regime where the 
distribution of fertilities [i is characterized by a power law ~ l//i 1+7 . For earthquakes, we expect 
such a power-distribution of fertilities with 7 = b/a based on the Gutenberg-Richter magnitude 
distribution ~ 10" bm and on the increase ~ 10 am of the number of aftershocks with the mainshock 
magnitude m. We derive the asymptotic distributions p r {r) and p g (g) of the total number r of 
offsprings and of the total number g of generations until extinction following a mainshock. In the 
regime 7 < 2 for which the distribution of fertilities has an infinite variance, we find p r (r) ~ l/r 1+ ~ 

and p g {g) ~ f- 1 . This should be compared with the distributions p r {r) ~ l/r 1+ 2 and 

p g {g) ~ l/g 1+1 obtained for standard branching processes with finite variance. These predictions 
are checked by numerical simulations. Our results apply directly to the ETAS model whose prefered 
values a — 0.8 and 6=1 puts it in the regime where the distribution of fertilities has an infinite 
variance. More generally, our results apply to any stochastic branching process with a power-law 
distribution of offsprings per mother. 

PACS numbers: 64.60.Ak; 02.50.Ey; 91.30.Dk 



INTRODUCTION 

All large earthquakes are followed by an increase of 
seismic activity known as "aftershocks." Aftershock se- 
quences of small earthquakes are less obvious because 
the aftershock productivity is weaker, but can be clearly 
observed when stacking many sequences 0. It is thus 
natural to assume that each earthquake can trigger its 
own aftershock sequence, and that observed aftershock 
sequences result from the cascade of direct aftershocks 
(triggered directly by the mainshock) and indirect after- 
shocks (triggered by a previous aftershock of the main- 
shock). This assumption is the basis of the Epidemic 
Type Aftershock Sequence model (ETAS) of seismicity 
[13 H EHH which describes earthquake trig- 
gering as a branching process. In addition, the ETAS 
model includes the Omori law decay ~ \/{t + c) p of the 
number of direct aftershocks as a function of the time 
t since the mainshock, where c is a small constant and 



p is an exponent close to but larger than 1. Previous 
works on this model have shown that the ETAS model 
provides a better fit to aftershock sequences than a sin- 
gle Omori law (no secondary aftershocks) || and that a 
significant fraction of aftershocks, of the order of 80%, 
are secondary aftershocks 0,0. The ETAS model has 
been used in many studies 0^0] to describe or pre- 
dict the spatio-temporal distribution of seismicity and 
reproduces many properties of real seismicity, including 
a renormalization of the Omori exponent from the lo- 
cal Omori law (direct aftershocks) to the global Omori 
law (observed rate of aftershocks including secondary 
aftershocks) [H |U HI , Bath's law a diffusion of 
aftershocks and realistic foreshock properties 0, Il3| . 
In this work, we present an analytical derivation of the 
distribution of the total number of aftershocks, summed 
over all generations of the cascade of aftershock trig- 
gering, and of the distribution of the total number of 
generations of aftershocks before extinction. 



There are two well known statistical laws that de- 
scribe the scale-invariance of earthquake physics with re- 
spect to magnitudes. First, the (complementary cumu- 
lative) Gutenberg-Richter (GR) distribution of earth- 
quake magnitudes gives the probability 

Pm (m) ~ 10- bm (1) 

that an earthquake has magnitude equal to or larger 
than m. This magnitude distribution p m (m) is not de- 
pendent on the magnitude of the triggering earthquake, 
i.e., a large earthquake can be triggered by a smaller one 

um 

Second, the average number of aftershocks triggered 
directly by an earthquake of magnitude m is found to 
increase with m as 

N m = K 10 am (2) 

where K is a numerical factor independent of the magni- 
tude. The number of direct aftershocks cannot be mea- 
sured, because what is observed is the total number of 
direct and secondary aftershocks. If earthquake trig- 
gering can be described by a branching process such as 
the ETAS model, then it can be shown that the scaling 
of the total number of aftershocks with the mainshock 
magnitude also obeys the law , but with a larger fac- 
tor K' which accounts for the cascades of secondary af- 
tershocks |lOj. The exponent a in can thus be mea- 
sured from a fit of the total number of aftershocks with 
the mainshock magnitude. Fits of the total number of 
aftershocks as a function of the mainshock magnitude in 
individual sequences su pport lt2Tl with an exponent a in 
the range 0.75-1 \A I HI l25l \?>\\ . However, the precision 
of these studies is limited by the narrow range of main- 
shock magnitudes considered and the large scatter of the 
number of aftershocks per mainshock. The value of a 
estimated in these studies may also be biased by the ar- 
bitrary constraint that aftershocks must be smaller than 
the mainshock, by the incompleteness of the catalog just 
after the mainshock, and by the background seismicity 
at large times after the mainshock. Other studies have 
measured the exponent a in 101 by calibrating the ETAS 
model to real data (individual aftershock sequences or 
complete catalogs) using maximum likelihood methods 
|M U-J Ell • These studies found a large scatter of the 
a exponent in the range 0.2 — 2. It is not yet clear if this 
range of values reflects a real variability of a or the inac- 
curacy of the estimation of a. One of us used a stacking 
method to estimate the average rate of earthquakes trig- 
gered (directly or indirectly) by a previous earthquake as 
a function of the magnitude of the triggering earthquake 
0, without constrain on the aftershock magnitude. For 
the catalog of the Southern California Data Center for 



Southern California, using the time period 1975-2003 
and m > 3 earthquakes, a is found equal to 0.8 ±0.1, 
smaller than b = 1 ± 0.1. Small earthquakes are thus 
collectively more important than larger earthquakes for 
earthquake triggering if a < b, because they are much 
more numerous than larger ones. 

Let us combine the two laws Q and |J2J to get the 
unconditional probability density for the number N m of 
events triggered directly by any event whose magnitude 
m is drawn at random from the GR law. For this, we 
note that 

Pr(fertility > N m ) = Pr(K 10 am > N m ) = Pr(m > - log 

a 

(3) 

where Pr means "Probability" . The first equality makes 
use of (0) and the third equivalence makes use of QJ. 
Hence, the probability density of the fertilities N m , for 
a magnitude m drawn at random in the GR law, is 

PN m (N m )~— , (4) 

with 7 = b/a. Typically, 6 = 1 and a w 0.8 leading 
to 7 w 1.25. Because 7 < 2, the variance of the fertility, 
for an earthquake of arbitrary magnitude drawn from 
the GR law, is mathematically infinite (if we assume the 
absence of a roll-off in the GR distribution, see below). 

Beyond earthquakes, a multitude of phenomena can 
also be described by branching processes with power-law 
distributions of fertility. Stochastic branching processes 
indeed describe well a multitude of phenomena [ij, |24j 
from chain reactions in nuclear and particle physics, 
material rupture, fragmentation and earthquake pro- 
cesses, critical percolation cluster sizes and population 
growth models, to population and biological dynam- 
ics, epidemics, economic and social cascades and so on. 
Branching processes are also of particular interest be- 
cause deep connections have been established with crit- 
ical phenomena |29l l30l | . Epidemic transmission of dis- 
eases, and more generally transmission processes involv- 
ing avalanches spreading on networks such as the World 
Wide Web, cellular metabolic network, ecological food 
webs, social networks, and so on exhibit such heavy-tail 
probability density functions (PDF) given by JSJ) below, 
as a consequence of the well-documented power law dis- 
tribution of connectivities among nodes. Our results 
are thus relevant to systems in which the number of 
offsprings may be large due to long-range interactions, 
long-memory effects or large deviation processes. Goh 
et al. actually derive results that overlap with ours 
in the context of avalanches in social networks. 

In branching processes with a finite variance of the 
number of daughters per mothers, various quantities ex- 
hibit power law distributions with universal exponents 
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at criticality (in statistical physics, the term "universal" 
refers to the independence of the critical exponents on 
the microscopic details of the physics) . This includes the 
distributions of cluster sizes, of the number of genera- 
tions before extinction and of durations which are mean 
field P, 01 ( "mean field" refers to the branching approx- 
imation which leads to a lack of dependence on the space 
dimension) . In the case of earthquakes and for other sys- 
tems mentioned above, the distribution of fertilities has 
an infinite variance, leading to an anomalous scaling of 
offsprings and generation numbers. While the number 
of direct aftershocks per mainshock for a fixed main- 
shock magnitude has a finite variance, usually modeled 
by a Poisson distribution, the effect of multiple cascades 
of triggering and the variability of the fertility of each 
earthquake lead to a much larger power-law distribution 
for the total number of aftershocks summed over all gen- 
erations. As a consequence, there are huge fluctuations 
of the total number of aftershocks from one sequence to 
another one, for the same mainshock magnitude. 

The goal of this paper is to provide a general ex- 
act derivation of the distribution of the total aftershock 
productivity of any earthquake, summed over all gen- 
erations of the cascade of aftershock triggering. We 
also derive the exact distribution of the total number 
of generations of aftershocks before extinction in the 
case of a power-law distribution of fertility relevant for 
earthquakes and for many other systems. Beyond earth- 
quakes our results apply to any branching process with 
a power-law distribution of fertilities. 



MODEL AND MAIN RESULTS 

We consider a general branching process in which 
each progenitor or mother is characterized by its av- 
erage number N m = Kfi(m) of children (first generation 
offsprings), where fi(m) = 10 Q (™ l_m °) is a mark asso- 
ciated with an earthquake of magnitude m > mo, k is 
a constant factor and tjiq is the minimum magnitude 
of earthquakes capable of triggering other earthquakes. 
We note /i the mark of an earthquake which has an ar- 
bitrary magnitude m drawn according to the GR law. 
According to @ , the mark /i is distributed according to 

7 

pM = 7TP7 , 1 < l" < 7 = b/a . (5) 

Note that is normalized: f^~°° d/i Pfi(/j) = 1. The 

relation N m = nfi(m) together with JSJ thus ensures 
that the fertility obeys the law (@J. For a fixed 7, the 
coefficient k then controls the value of the average num- 



ber n of children of first generation per mother: 

7 

n = (N m ) = k{/j,) = k , (6) 

7 - f 

where the average (N m ) is taken over all mothers' mag- 
nitudes drawn in the GR law. In the terminology of 
branching processes, n is called the branching ratio. For 
n < 1, there are less than one child per mother: this 
corresponds to transient (sub-critical) branching pro- 
cesses with finite lifetimes with probability one. For 
n > 1, there are more than one child per mother: this 
corresponds to explosive (super-critical) branching pro- 
cesses with a number of events growing exponentially 
with time. The value n = 1 of exactly one child per 
mother on average is the critical point separating the 
two regimes. 

The realized number of children of an earthquake of 
fixed magnitude m can be deterministic or may result 
from a Poisson or other more general distribution with 
mean Kfi(m). Any given child i may then generate an 
average number n/Ji of children, where the mark /i^ is 
specific to the child i and is drawn from the PDF JSJ. 
These grand-children of the initial progenitor in turn 
generate new children with marks drawn from JSJ, and 
so on. The process cascades down along generations. 

Here, we focus on global quantities p r (r) (total num- 
ber r of offsprings until extinction following a main- 
shock) and p g {g) (total number g of generations un- 
til extinction following a mainshock). Therefore, arbi- 
trary time-dependent branching processes can be con- 
sidered and our results thus apply to arbitrary stochastic 
marked point-processes in discrete or continuous time. 
In particular, our results apply to the ETAS model. We 
have previously shown that, for n < 1, the infinite vari- 
ance of the number of first-generation events leads to 
anomalous global direct and inverse Omori law with ap- 
parent exponents varying continuously with 7 .13]. Our 
results also apply to variations of the ETAS models with 
different arbitrary time dependences. 

In this work, we present an analytical derivation of the 
distribution of the total number of aftershocks, summed 
over all generations of the cascade of aftershock trigger- 
ing, and of the distribution of the total number of gener- 
ations of aftershocks before extinction. Specifically, we 
uncover a novel regime with continuously varying expo- 
nents for the probability density function (PDF) p r (r) of 
the total number r of progenies and p g (g) of their total 
number g of generations of aftershocks of a mainshock 
before extinction. This regime appears in the regime of 
large deviations 7 < 2 relevant for earthquakes, when 
the distribution of the number of first-generation off- 
springs from any mother has a power law tail with infi- 
nite variance. 
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We study the sub-critical and critical regimes for 
which the number n of children per mother, averaged 
over all possible numbers of children per mother, is less 
than or equal to 1. This condition n < 1 ensures that, 
with probability 1, the cascade ends after a finite num- 
ber of generations with a finite total number of offsprings 
[lll24j. Figure ^ presents comparisons between numeri- 
cal simulations of branching processes for different val- 
ues of 7 and our main result (at criticality n = 1) derived 
below: 

p r {r) « C r /r 1+ ^ ; Pg (g) « C g /g 1+ ^ , (7) 

for 1 < 7 < 2. C r and C g are two constants independent 
of r and g, respectively, such that p r (r) and p g {g) are 
normalized to 1. For n < 1, the intermediate asymp- 
totics @ holds for r < R 2 ^ 1/(1 - n)^" 1 ) beyond 
which there is another power-law asymptotic 

p r (r) ~ l/r 1+ ^ . (8) 

For 7 > 2 and n = 1 , the standard mean field scaling 

Pr (r)~l/r 3 / 2 and p g (g) ~ l/ 5 2 (9) 

arc retrieved. The regime < 7 < 1 gives rise to explo- 
sions and, in continuous time, to stochastic finite-time 
singularities [2^] . If the GR law is truncated or exhibits 
a deviation from its standard form in the large magni- 
tude range, our results will hold for intermediate values 
of r and g and will be truncated at large r's and large 
g's. 

The following section gives the technical derivation of 
the results {ZEE} . 

GENERAL FORMULATION 

Formal solution in terms of Probability Generating 
Functions (PGF) 

The most general relations of branching theory are ex- 
pressed via probability generating functions (PGF) de- 
fined by 

00 

(* fl )=5>r(r)* p (10) 

r=0 

where (...) means statistical averaging. By definition, 
p r (r) — Pr(i? = r) is the probability that the random 
variable R takes the value r. It can be obtained from 
its PGF by the relation 

1 d r (z R ) 

p r (r)=Pr(R = r) = --±-<- . (11) 

T. CLZ ~ n 
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FIG. 1: Numerical tests of the anomalous scalirigs Q at 
criticality n = k{(j) = 1. The noisy lines are the Monte- 
Carlo simulations. The smoothed lines are the predictions 
Q. The deviation between theory and numerical simula- 
tions on P g (g) for 7 = 1.1 results from the very large theoret- 
ical exponent equal to 1 + — Lj- =11, describing a fast fall-off 
which is hard to document numerically since one decade in 
the horizontal scale corresponds to 11 decades on the vertical 
scale. 

The key property of PGFs is that the PGF of the sum of 
statistically independent summands is equal to the prod- 
uct of the summands PGF's. This is useful for branching 
processes for which the number of daughters triggered 
by different mothers are statistically independent. We 
introduce four PGF's: (i) G^ (z) is the PGF of the num- 
ber i?i of daughters generated from a given mother with 
fertility k/i at the first generation; (ii) G(z) is the same 
as Gjj,(z) but averaged over all mothers' fertility n/i; (iii) 
M (z) is the PGF of the total number of daughters gen- 
erated from a given mother with fertility Kfi summed 
over all generations; (iv) Q(z ) is the same as M (z) but 
for a mother of arbitrary fertility. Then, general branch- 
ing theory 0, [24| gives the functional equations 

e(z) = G(z&(z)) ; 9p[z)=G ll {zQ{z)). (12) 

Using (R) = d@{z) / dz\ z= \ and similarly for the other 
PGF, the formulas (JEJ) lead to (R) = n/(l - n) for the 
number of daughters summed over all generations and 
averaged over all possible mother's fertilities, where n 
defined by 10 is the average number of first-generation 
daughters from mothers of arbitrary fertility. 

Using Lagrange series, we transform the implicit equa- 
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tion ((12(1 in M into an explicit equation. Recall that a 
Lagrange series is the Taylor expansion of the function 
Gn(y(z,a)) with respect to z, where y{z,a) is solution 
of the implicit equation y — a + zG(y). For infinitely 
differentiable functions G M (y) and G(y) in the neigh- 
borhood of z — a, we have the following Taylor series 



G,(y) = G,(a) + J2^ 



fc=i 



kl da*- 1 



G\a) 



da 



(13) 



In H12|) . we have y(z,a — 0) = z0(z), and using the 
identity p ri/i (r) = (l/r!)cT6 A1 (z)/dz r | z= o allowing to re- 
cover the PDF from its generating function, we obtain 
the explicit formula for the PDF p r ,^(r) of the total 
number of daughters born from a mother with fertility 



Pr,„(0) = G M (0) 

i dr- 

P ^ {r) = Til? 



G r {z) 



dGjjjz) 
dz 



(14) 
r >(Q5) 



z=0 



Let us now specialize to the Poisson distribution 



(/^) ri 



(16) 



giving the PDF of the number i?i of daughters of first 
generation born from a mother with fertility fi. By con- 
struction, the average of R\ at fixed mother fertility /j, 
over an ensemble of Poisson realizations is = l^n. 

The PGF associated with JTJJ is 



^k(z-1) 



(17) 



From 1(17(1 . we obtain the expression of G(z) by av- 
eraging over all possible fertilities fi according to the 
PDF JSJ) . Using explicitly the normalized PDF (fi) = 
7//X 1 " 1 " 7 for (ti > 1 and otherwise, we obtain 



G{z) = 7 «T(1 - z) 7 r(- 7 , k(1 - 2)) . 



(18) 



Expression (|18fl can be expanded as 

G(z) ~ 1 - ai y - a 7 y 7 + a 2 y 2 + ... (19) 



where y = k(1 — z), a\ = t^j, a 7 = T(l — 7) and 
a 2 = 2 {-y-2) ■ This expansion $B^ , valid for any 7, ap- 
plies beyond the specific Poisson process ((16(1 • It is solely 
based on the asymptotic power law PDF (JjJJ of the aver- 
age number of daughters at the first generation. Thus, 
our asymptotic results Q derived below hold under very 
general conditions. 



Probabilistic interpretation 

Let us provide a probabilistic interpretation of for- 
mula HI 5(1 . In addition to be of intuitive appeal, it will 
be very useful for the derivation of our results. 

Let us consider the random integer R with PDF 
P(k) = Pr(R = k) and its PGF as 



5(z)^^P(fc)z fc , 



k=Q 



and reciprocally 



P{k) = Pr(JZ = k) = 1 



k\ dz k 



(20) 



(21) 



2=0 



Let us decompose R = X + Y, where X and Y are sta- 
tistically independent random integers with PGFs Q{z) 
and G(z) respectively. Due to the statistical indepen- 
dence of X and Y, we have S(z) = Q(z)G(z) and rela- 
tion ((21(1 takes the form 



1 d h 



k\ dz k 



: [Qiz)G{z)\ 



P{k) = Pr(X+Y = k) 
Let in turn 

Y = Y{r) = Yy + Y 2 + ■ ■ ■ + Y r 



z=0 



(22) 



(23) 



where {Yi,.. . ,Y r } are r statistically independent ran- 
dom integers with the same PGF G(z). Then, formula 
J23 takes the form 



P(k) = Pr(X + Y(r) = k) = I A_ [Q(z)G r (z)} 



z=0 

(24) 



Let us now rewrite relation 1(151) in a form similar to 
1(24(1 . For this, let us introduce the auxiliary PGF 



dG^z) 



(ri(n)) dz 



(25) 



It is easy to prove rigorously that, for an arbitrary PGF 
G^{z) possessing a finite mean value 



dGfi(z) 



dz 



< 00 



(26) 



z=l 



then the auxiliary function Q M (z) is indeed the PGF 
of some random variable. Moreover, in the case under 
consideration of a Poissonian PGF, we have 



G^z) = e^ (z - 1} 
and 



) fjt ( z ) = G IM (z)=e^-^ 
(27) 



(28) 
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One can thus rewrite relation l|15|) in the form 



(29) 



z=0 



Taking into account relation (|24|l and introducing the 
new random integer X(kh) possessing the Poissonian 
PGF (J23, we can rewrite (J23 in the form 



PrAr) = — Pr [ X M + Y ( r ) =r-l] 



(30) 



This is the key formula that we will use for deriving 
the asymptotic relations for the PDF p r ll (r) using limit 
theorems. 



DISTRIBUTION OF THE TOTAL NUMBER OF 
AFTERSHOCKS 



Case 7 — » oo 

Let us first consider the limiting case 7—5-00. The 
formal limit of a power law with an exponent going 
to infinity can be assimilated to an exponential in the 
following sense. Writing the tail of the complemen- 
tary distribution as A 1 /(A + ^) 7 , and assuming that 
A also grows with 7 as A = 7/a, then A 1 /{A + ^) 7 = 
1/(1 + a^/7) 7 ^ 7 ^+oo exp[— afi\. This corresponds to 
an exponential tail for the distribution of the number of 
first generation daughters for arbitrary mother's fertili- 
ties, as in l)16f) but with the same fertility for all moth- 
ers. In this case, n = n as can be seen from This 
amounts to replacing expression (|18(l by G(z) — e K ( z-1 ' 
and we also have G fl (z) = e M "( 2_1 ). Thus, 

G r (z) ^(f) =/iKe «(*-l)(H-M) . 



dz 



Correspondingly, 



dz r 



G r (z) 



dG^z) 
dz 



(31) 

By substituting 1)31(1 into the right-hand-side of (|15[1 
gives exactly 



fV( r ) = ^n r e 



- Kfi \(P + r)e K ] 
r\(n + r) 



(32) 



There are two ways of deriving asymptotic formulas 
for p r A r ) m the case 7 = 00. The first one relies on 
Stirling's formula r! ~ \j2itr r r e~ r applied to expression 
(|32|) . The second method relies on the probabilistic in- 
terpretation of the general formula (|15|) given in section 
□ together with the central limit theorem. The formulas 



obtained by these two methods differ in details, but are 
actually very close quantitatively over the most part of 
the tail, and are easily checked to converge to the same 
result asymptotically. The approach using the Stirling 
formula works only for 7 = 00 and does not allow to 
derive more general results for 1 < 7 < 00. The sec- 
ond approach is much more powerful and elegant. In 
particular, all asymptotics and crossovers in the case 
1 < 7 < 00 are obtained by the second approach. 

Let us first examine the first method. Using the above 
mentioned Stirling's formula, one can rewrite 1)32(1 in the 
form 



1 r 
PrA r ) = ~ f( x ) x = - 



where 



c 



'2ttx (1 + x) \ x 



fix 



and 



C = yfjxe 



Kfl 



In/ 



1 



(33) 



(34) 



(35) 



Because Stirling's formula is very precise even for r ~ 1 
(for example 1! = 1 while Stirling's formula gives l! ~ 
0.922, 2! = 2 while Stirling's gives 2! ~ 1.919 and so 
on), formula ((34(1 actually works well even for not too 
large /x and intermediate n < 1. The advantage of ((34(1 
is that it gives a clear understanding of the structure of 
the PDF f(x). It consists of (i) a characteristic power 
tail 



r-3/2 



(ii) an exponential decaying factor 
f(x) oc e~» Sx 



(36) 



(37) 



which disappears when k = 1, i.e., 5 = 0, and (iii) an 
algebraic factor 



1 \ f^ x 

il + -J -'"cxp( ^) (38) 



which possesses a lower cut-off at x < /.i. 

The shortcoming of l(34|) is that we can derive it only 
when an explicit expression such as ((32(1 is available. 
We have not such luxury in the more general case 1 < 
7 < 00 for which we need the second approach. Let us 
first quote the asymptotic formula ((39(1 given below for 
7 — * +00 obtained from the application of the second 
probabilistic method. We shall show how to derive ((39|> 
as a special case of the next section for 2 < 7. 
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When the variance a\ of the number r% of daughters 
of the first generation from mothers of arbitrary fertili- 
ties is much greater than 1 — re, we find (see next section) 
that, due to central limit theorem, expression (|32|) re- 
duces to 



Pju( r ) 



fjLK 



r^/2n re (r - 
At criticality, n = 



■A*) 



exp 



[(1 



jJLK 



2re (r + fi) 

(39) 

this expression becomes 
exp (-£) , (40) 



which retrieves the announced well-known mean field 
asymptotics p^(r) ~ 1/r 3 / 2 . The exponential term 
exp(— /i 2 /(2r)) describes the roll-off of the number of 
aftershocks for small numbers close to the characteristic 
value r ~ fx . 

In the subcritical regime n < 1, expression Q32I ) gives 
an exponential decay of the number of aftershocks for 
large r. 

These results can be checked explicitly as follows. The 
case of 7 — > +oo corresponds either to b — > oo (all events 
have the same magnitude) or to a — (all the events 
trigger the same expected number of offsprings indepen- 
dently of their magnitude). In both cases, \x tends to a 
constant independent of the mainshock magnitude. In 
order to check the above results, let us take this con- 
stant equal to 1. This choice will modify the constants 
but not the functional dependences. Assuming [i = 1 
transforms ilML'l) into Taking jj, = 1 transforms l|32|) into 



Pfj.ir) = re' 



(r + 1)' 



-«(r+l) 



(41) 



(r + 1)! 

By using the Stirling formula, it has the approximation 



Q r In K.+ (r+l)(l — k.) 

\/2^ (r + l) 3 / 2 



(42) 



When re = n = 1 (critical case), we retrieve PuO") ~ 
r -3/2_ when re = n < 1 is close to 1, p^(r) ~ r~ 3 / 2 e~ 5 r 



where <5 is defined in (|35|) . 



Case of finite variance 2 < 7 < 00 

In this case, the summands in the sum <|23|) have a 
finite mean (Y{) and finite variance a\. Thus, the Cen- 
tral Limit Theorem (CLT) holds (see for instance (2?| 
for a pedagogical exposition of the CLT). Let us recall 
the explicit expressions of the mean and variance: 



(Yx) =n = 



re7 



7(7-2) 



+ n , 



(43) 



If r ^> 1 (in practice, it is sufficient that r > 6), the sum 
Y(r) in (|23|l converges to a Gaussian variable with the 
following mean and variance 



(Y(r)) = r n . 



cr 2 (r) = ra\ 



(44) 



Let /ire be integer (just for illustrative purpose). Then, 
one can imagine X(//re) as a sum 



A(/ire) 



f I hi 

m— 1 



of independent summands {A"i,A"2, 
following PGF, mean and variance: 



(45) 



, X^ K } with the 



Q{z) 



(Xi) 



(46) 



If the number of summands /ire 3> 1, then the sum (|45l) 
is asymptotically Gaussian as well, with mean and vari- 
ance equal to /ire. 

Thus if r 3> 1 and /ire ^> 1, then the sum X(/j,K)+Y(r) 
is asymptotically Gaussian and we can use the following 
asymptotic formula 



Pr [A(re/i) + Y(r) 



exp 



(m - a M (r)) 2 
2a 2 (r) 

(47) 

Here (r) and cr 2 (r) are respectively the mean and vari- 
ance of the sum A(/ire) + Y(r): 

an(r) = rn + fiK , a M (r) = + //re . (48) 

Substituting (gSJ into (HTJl and E7|) into l{3TJj) . we obtain 



Pr^O") 



//re 



T\Jl'n(a\r + /ire) 



exp 



[(1 — n)r — /ire — 1] 
2(er 2 r + /ire) 



(49) 

For 7 = 00, for which we have the mean and vari- 
ance defined in (|43[) reduce to a\ = n = re, relation 1(49(1 
transforms into the previously announced result l|39|) . In 
this case, we can test the quality of the general asymp- 
totic relation (|49J) by comparing its particular applica- 
tion 1(39(1 with the exact expression 1(32(1 and its Stirling 
asymptotics (|34|l . The exact and Stirling's approxima- 
tion essentially coincide while the "Gaussian" approxi- 
mation is also excellent and goes closer and closer to the 
exact formula the larger // is and the closer n — re is to 
1. 

For o\r > /ire, expression 1)49(1 can be simplified into 



PmO) 



//re 



exp 



[(1 — n)r — /Ltre — iy 
2afr 



(50) 
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which reduces again to the standard mean field asymp- 
totics (|40|l at the critical point n = 1 and for r 3> [i 2 K 2 . 
Expression (|50|) is obtained using the following approx- 
imation 



[IK 



a 2 r . 



(51) 



that is, by neglecting the term fin. From a probabilistic 
point of view, this corresponds to using the Law of Large 
Numbers for the sum (|45|l . Namely, if /xk 3> 1, then one 
can replace (14511 by its mean field limit 



X([lk) ~ (X([ik)) = [IK 



(52) 



which provides the truncated equality (|51() . Another 
reason for using the mean field limit ((52(1 is more tech- 
nical: in order to be able to effectively use the same 
technique in the case 1 < 7 < 2, we need to implement 
such a "mean field approximation." Otherwise, the for- 
mulas would include very complicated convolutions of 
Gaussian and Levy stable laws. 

We have checked the validity of the mean field approx- 
imation by comparing numerically the general asymp- 
totics ((49(1 and of its truncated mean field limit version 
(|5Tj|) for various values of 7, n and [ik. We find always 
an excellent convergence in the tail. For instance, con- 
sider the simplest case 7 = 00, for which the general 
asymptotic formula 1)49(1 transforms into 1)39(1 . while its 
truncated mean field version is 



PrA r ) 



[IK 



exp 



[(1 - k)t - [ik- If 
2nr 



(53) 

The agreement between all these expressions is excel- 
lent in the tails for all values of the parameters. The 
truncated mean field approximation describes satisfac- 
tory the body of distribution and becomes very precise 
in the tail of the distribution, even for moderate values 
of /i 25 and n = n < 0.9. Even better agreement 
is obtained for larger [i and closer to the critical case 
n = 1. 

Let us complement this section by a few additional 
useful formulas. The expansion ((19(1 with = 
dG(z) I dz\ z= i shows that the average number r\ of first- 
generation daughters from mothers of arbitrary fertili- 
ties is n = ([j)k = kj/(j — 1) and that, for 1 < 7 < 2, 
the variance of r\ is infinite. For 7 > 2, the variance of 
r\ is o\ = n + n 2 /[j(j — 2)]. The form of 1(19(1 is essen- 
tially controlled by the expression of the characteristic 
function of the PDF JSJ) and it is thus not a surprise 
that the PDF of the number of daughters of first gener- 
ation born from a mother with arbitrary fertility has the 
same asymptotic form as (0. For a mother with fixed 
fertility /it, the average and variance of the total number 



of daughters are respectively 
[XK 2 M K 



1-n ' " (1-n) 3 



7(7-2) 



for 7 > 2 . 
(54) 



Case of infinite variance 1 < 7 < 2 



We now turn to the novel regime 1 < 7 < 2. In this 
case, it is convenient to rewrite relation 1(15(1 in a more 
transparent probabilistic form 

r 

P„(r) - — Pr [X + Y(r) = r - 1] , Y(r) = ^Y fe , 

k=l 

(55) 

where {X, Y\, . . . , Y r } are mutually independent random 
integers such that the PGF of X is G M (z) given by ((1711 
while the PGF's of the remaining integers are equal to 
G(z) given by 1(18(1 . For 1 < 7 < 2, the variance of 
each variable Yfc is infinite, for the same reason that 
a\ is infinite. From the generalized limit theorem [l8|. 
Y(r) — Y(r) — nr converges in distribution to a PDF 
f y (s) which is proportional to a stable infinitely-divisible 
PDF (f-y(x) with exponent 7 



fv(s) 



-1/7 



¥>7 



-1/7 



where 



= «br(-7)] 



1/7 



(56) 



(57) 



and 



¥>7 



OO 

(x) ~ J exp (^u 1 cos (~ cos sin(— -^-) + u x^j du 



(58) 

with <p 7 (0) = — l/r(— 1/7). Thus we obtain the distri- 
bution of aftershock number from ((55() and l(56() 



[IK ( (1 — n)r — uk — 1 
Mr) - ~T+T ^7 -177 



(59) 



For large number of aftershocks r >> Ri where R\ is 
defined by 



Ri 



( [IK + 1 



(60) 



\ v ) 7IX-7) ' 
we can use the following asymptotic of Lp 1 (x) for x — > 00 



<P-y(x) 



-1-7 



r(-7) 



(61) 
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10" 10' 10" 
total number of events 



FIG. 2: Numerical tests of the cross-over between Q and 
llrjit which is predicted by lEU of the PDF p M (r). The thick 
(resp. thin) line is the numerical simulation (resp. prediction 
H59^ ) for n = 0.97 and 7 = 1.5. The two dashed lines outline 
the first intermediate asymptotic with exponent 1 + I/7 for 
-Ri -C r >C R2 and truly asymptotic with exponent 1 + 7 for 
r ^> i?2- The crossovers Ri — 1.5 and R2 = 1.6 x 10 4 are 
defined by (I60H and 16311 respectively. 



Putting (|61fl in (|59|) retrieves at criticality (n = 
1) our main result announced in Q with C r — 

/ i^(o)/[ 7 r(- 7 )] 1 /7. 

In the subcritical case (n < 1), there is another power 
asymptotics. For r R\, one can rewrite l|59|) in the 
form 



Pn(r) 



1 



n 7-i 

T T 



[ 7 r(-7)]V7 r /3 • v v 

If additionally r ^> R 2 , where R 2 is defined by 

7/(7-1) 

R2 = ' 



(62) 



1 



(63) 



then, using the asymptotic (|61|l of ip-y(x), we obtain 



Pn(r) — A* 7 



"(7- 1) \ 
(l-n)7/ 



7+1 



,1+7 



(64) 



for r ^ max{i?i, R 2 }. If additionally i?2 3> -Ri, which 
holds if ^ < (7 - l)r(- 7 ) ^-t, then for i?i < r < 
i?2j there is an intermediate asymptotics following ex- 
pression Q). These results are checked in Figure |21 



DISTRIBUTION OF THE TOTAL NUMBER OF 
GENERATIONS 

We now turn to the determination of the PDF of the 
total number of branching generations, in other words, 



to the probability that the branching process terminates 
at a given generation number. Let U^{z; g) (respectively 
U(z\g)) denote the PGF corresponding to the number 
of daughters born from a mother with fertility /i (respec- 
tively with arbitrary fertility) at the 5-th generation. A 
standard result of branching theory is |l| 

U fl (z;g + l) = G^(U(z;g)) ; U(z;g + 1) = G(U(z;g)), 

(65) 

where U(z; 0) = z. The probability that the branching 
process survives at the g-th generation is 



JV(S) = 1-^(0, g) 



(66) 



The probability of termination of the branching process 
at the g-th generation is then given by 



p a Aa) =pM -Puig + i) ■ 



(67) 



From H65I66JI . the probability for the branching process 
to survive at the g-th. generation is 



p^g + l) = l-G^l-p{g)). 



(68) 



p(g) is the surviving probability at the g-th generation 
triggered by a mother of arbitrary fertility //, which 
obeys the recurrence 



p(g + l) = l-G(l-p(g)). 



p(0) = 1 



(69) 



Equations l|68[) and l|69(l are easily solved numerically. 
We now extract the asymptotic power law behavior of 
Pg,tj,{g)- For small p(g), the right-hand-side of (|69|l can 
be expanded as 



1 — G(l — p) ~ np - 



C 
7-1 



pi - D p 2 , (70) 



where C = rf> T(2 - 7) [(7 - l)/^ 7 and D = n 2 (7 - 
l) 2 /[27(7 - 2)]. For 7 > 2, it is enough to take the 
leading behavior of (|70|l up to second order in powers 
of p: 1 - G(l - p) ~ np - Dp 2 . With Jgj), this gives 
p(g + l)— p(g) = — (1 — n)p(g) — Dp 2 . Close to criticality 
n — ► 1 for which a large number g of generations occur, 
the leading behavior of the survival probability can be 
obtained by taking a continuous approximation to 169f) . 
giving dp ^ — — (1 — n)p — D p 2 . Its solution with initial 
condition p(go) = po has the form 



p(g) = 



p (l - ri) 



(1 



p D)e 



(l-")(s-So) 



PoD 



(71) 



Expression l|71|l gives a power law p{g) ~ t/{D g) for 
1 <C g < 1/(1 — ri) and crosses over to an exponential 
law for g > 1/(1 — ri). Knowing p(g) , p^ (g) given by (jrjS|) 
can be obtained. For instance, in the case of a Poisson 
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distribution (|16fl giving l|17|). we obtain p fl (g + 1) = 1 — 
exp [— /z np(g)]. For [inp(g) <C 1, P M (.9 + 1) becomes 
£V(g + 1) ~ nnp{g). The sought PDF of p g ^(g) defined 
by (|u7|l is given by 



PgAd) - M« 



dp{g) 

dg ~ D g 2 



(72) 



The last equality holds for g > 1/(1 — n). This retrieves 
the standard mean field asymptotics p g ,^,{g) ~ l/ff 2 - 

For 1 < 7 < 2, it is sufficient to keep only the following 
terms 1 — G(l —p)~np— -^jP 1 in the expansion lf7U|) . 
Taking again a continuous approximation of Ij69|) gives 
= —(1 — n)p — ^3jp 7 , whose solution, for g <C 



dg 

l/[(7 



1)(1 



l- 



n)], reads 
+ C( 5 



fhh 



-1/(7-1) 



c 



1/(7-1) 



(.9 - 9o ) 



(73) 

Correspondingly, the PDF Pg^ig) that the number of 
generations is exactly g is given by Q with C g = 
Hi (7-1/(7-1) ) va iid f or 7 not too c i ose to 2 . 



CONCLUDING REMARKS 

We have shown that the existence of cascades of trig- 
gered seismicity produces huge fluctuations of the to- 
tal number of aftershocks and of the total number of 
generations from one sequence to another one for the 
same mainshock magnitude, characterized by power-law 
asymptotics. If the distribution of offsprings per moth- 
ers (of direct aftershocks per mainshock for seismicity) 
has a finite variance, we recover the well known mean- 
field results (2J. In the case of an infinite variance rele- 
vant for earthquakes, we have discovered a new regime 
with exponents that varies continuously with the expo- 
nent 7 = b/a. The anomalous scaling reported here for 
1 < 7 < 2 gives rise to less wild fluctuations in the to- 
tal number of daughters from one mother to the next, 
compared with the mean field regime 7 > 2. For in- 
stance, for earthquakes, we have probably 7 sa 1.25 for 
the preferred values b ss 1 and a ss 0.8 [9J, which leads 
to p r (r) ~ i/yi+0-8 anc j p g (g} ^ 1 / g 1+i at 7i = 1 , com- 
pared with p r (r) ~ l/yi+O- 5 and p g {g) ~ l/g 1+1 for 
7 > 2. The reason for this behavior lies in the impor- 
tant role exerted by the mothers with the largest fertil- 
ities oc fx on the rate of daughter births. In particular, 
this role explains why, away from criticality but close 
to it (n < 1), p r (r) crosses over from l/r 1+ ~ to l/r 1+7 
as r increases: the latter asymptotic regime is nothing 
but the tail of the PDF (JSJ) which controls the number 
of daughters born in the first generation from a mother 
with arbitrary fertility. 



How relevant is the infinite variance regime 7 < 2 in 
view of the uncertainties in the exponents a and &? In- 
deed, the value of a is not only approximately known 
but its constancy (as a function of time) and uniformity 
(as a function of space) remains to be ascertained. In 
addition, the 6-value also fluctuates and is typically in 
the range 0.75 — 1.25 Q. If we take the largest &-value 
of this interval (b = 1.25), we remain in the infinite vari- 
ance regime 7 < 2 as long as a > 0.63. Most present 
estimations discussed in section 2 give values above this 
lower bound. It thus seems that the results presented 
here should remain pertinent as more precise estima- 
tions of a and b become available. 

All these discussions assume untruncated power law 
distributions. However, it is well-known that the 
Guten berg - Richter law exhibits an upper magnitude 
cut-off 0,H3] ■ By the logic leading to jtj, this automat- 
ically implies also an upper fertility cut-off. This in turn 
removes the divergence of the variance for 7 < 2. How 
does this affect the distributions p r (r) and p g {g)l This 
question is standard in the theory of power laws (see for 
instance H?} and references therein). Since the upper 
cutoff is very large so that the actual range of events 
is very large, our results obtained without truncation 
hold for a large range of values of r and of g respec- 
tively. However, the predicted power laws (Q have to 
cross-over to the mean-field ones (0 at threshold values 
beyond which the finiteness of the variance of fertilities 
become flagrant. In contrast, the presence of a lower 
cutoff magnitude 0] does not modify our results on the 
power law tail and has an influence only in shaping the 
bulk (small values) of the distributions. 

The direct validation of our predictions Q and (JHJ 
on earthquake catalogs is not feasible, due to the im- 
possibility to distinguish aftershocks from uncorrelated 
events at large times after the mainshock, and due to 
the limited number of large aftershock sequences. Our 
results have however some consequences for the statisti- 
cal properties of aftershock sequences. We have indeed 
shown in that the existence of large fluctuations 
of the number of aftershocks per mainshock (summed 
over all generations) induces a non-trivial scaling of the 
difference in magnitude between a mainshock and its 
largest aftershock, so that Bath law can be recovered in 
the regime b > a. These results also suggest that the 
common use of a Poisson distribution in seismicity fore- 
casts to model the distribution of the number of events 
within a finite space-time window is questionable, since 
the simple physics of cascades of earthquake triggering 
gives a very different (power law) distribution. Recent 
observations show non-Poisson power law distributions 
of seismic rates (see [z| and work in progress). We sug- 
gest that these observations and our results could be 
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used to improve earthquake forecasting by providing a 
more realistic distribution of the number of events. 
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